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Abstract 

The nitrogen-vacancy center in diamond is a promising candidate for realizing the spin qubits 
concept in quantum information. Even though this defect is known for a long time, its electronic 
structure and other properties have not yet been explored in detail. We study the properties of the 
nitrogen-vacancy center in diamond through density functional theory within the local spin density 
approximation, using supercell calculations. While this theory is strictly applicable for ground state 
properties, we are able to give an estimate for the energy sequence of the excited states of this 
defect. We also calculate the hyperfine tensors in the ground state. The results clearly show 
that: (i) the spin density and the appropriate hyperfine constants are spread along a plane and 
unevenly distributed around the core of the defect; (ii) the measurable hyperfine constants can be 
found within about 7 A from the vacancy site. These results have important implications on the 
decoherence of the electron spin which is crucial in realizing the spin qubits in diamond. 

PACS numbers: 71.15.Mb, 71.55.Ht, 61.72.Bb, 61.72.Ji 
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I. INTRODUCTION 



The nitrogen-vacancy (NV) center in diamond has attracted much attention in recent 
years, because it has been shown to give rise to a single optically active level within the 
diamond band gapi^ , and as such provides an interesting candidate for a qubit for quan- 
tum computing applications.- 1 ^^^ Besides providing a single photon source for quantum 
cryptography,-*^ the NV center is also a promising candidate as an optically coupled quan- 
tum register for scalable quantum information processing, such as quantum communication^ 
and distributed quantum computation.— In addition, it has been recently demonstrated that 
proximal nuclear spins can be coherently controlled via hyperfine interaction^ and used as 
a basis for quantum memory with an extremely long coherence time.— Therefore, knowing 
the electron-nucleus hyperfine interaction and its position-dependence is essential to analyze 
and optimize coherent control of proximal nuclear spins.— 

Experimentally, the hyperfine constants of the closest atoms near the vacancy are known 
from electron paramagnetic resonance (EPR) and electron-nuclear double resonance (EN- 
DOR) studies.— 1 ^ The hyperfine interaction between 13 C isotopes farther from the vacancy 
contributes to the coherent electron- nuclear spin states in the measurements.— 1 ^ Ah initio 
supercell calculations can be a very useful tool for determining the hyperfine tensors of a 
defect. For instance, such calculations have been used to identify the basic vacancy defects 
in silicon carbide by comparing the measured and calculated hyperfine constants.— ii&22i2i In 
the present paper we focus on the calculation of the full hyperfine tensor of the NV center 
in diamond which is of very high importance for qubit applications. Previous theoretical 
work has reported a calculation of the hyperfine constants of the NV center in a small 64- 
atom supercell,— but that work determined only the Fermi-contact term rather than the full 
hyperfine tensor. We will show in the Results and Discussion section that the conclusions 
based on that earlier analysis were adversely affected by the small unit cell size. The larger 
supercell of 512 atoms employed here does not suffer from this limitation and provides a re- 
alistic picture for the defect properties, which are in excellent agreement with experimental 
measurements. 

The rest of this paper is organized as follows: Section II gives a general discussion of 
the electronic states of the NV center based on a single-particle picture and the many-body 
states that can be constructed from this basis. Section III describes the method of the 
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first-principles calculations we performed. Section IV presents and discusses our results 
concerning the atomic and electronic structure of the defect, as obtained from the ab initio 
calculations. In Section V we present a detailed discussion of the calculated hyperfine 
interactions. Finally, we give our conclusions on the nature of this defect in Section VI. 

II. THE ELECTRONIC STATES OF THE NV CENTER IN DIAMOND 

The electronic structure of NV center in diamond has been discussed in detail in a recent 
paper 23 ; we briefly review the main points here. The NV center was found many years ago in 
diamond. 24 The concentration of NV centers can be enhanced in N-contaminated diamond 
by irradiation and annealing.— >2£ The model of the NV center consists of a substitutional 
nitrogen atom adjacent to a vacancy in diamond.— i2ii2£i2li The NV center has a strong optical 
transition with a zero phonon line (ZPL) at 1.945 eV (637 nm) accompanied by a vibronic 
band at higher energy in absorption and lower energy in emission. Detailed analysis of the 
ZPL revealed that the center has trigonal, C3„ symmetry.— Later, an optically induced EPR 
center was found in diamond which correlated with the NV center.™ The EPR center showed 
trigonal symmetry with a spin polarized triplet state (S=l). Since the nitrogen atom has 
five valence electrons and the S=l state implies even number of electrons, the NV defect 
must be charged in the EPR measurement. It was assumed that the NV defect is negatively 
charged and the extra electron may be donated from isolated substitional nitrogen defects.— 
In a recent measurement, the coupling between the NV center and the nitrogen substitional 
has been indeed detected.— Loubser and van Wyk— measured the NV EPR signal just 
under the optical excitation, based on which they proposed that the spin polarization arises 
from a singlet electronic system with inter-system crossing to a spin level of a metastable 
triplet. Redman and co-workers™ detected the NV center in the dark even at 100 K by EPR, 
from which they concluded that the S=l state is the ground state of the NV center. Hole 
burning,™ optically detected magnetic resonance 30 and Raman heterodyne measurements 3 - 1 " 
also showed the S=l state to be the ground state of the NV center. 

A group theory analysis based on a single-particle picture can be very useful in under- 
standing the nature of the defect states and the possible optical transitions between them. 
While the number of electrons in the NV center has been disputed in the literature^ 3 - 3 - 
a previous ab initio calculation clearly supported the negatively charged NV defect^ 4 - as 
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was originally proposed by Loubser and van Wyk. 16 We will also show in the Results and 
Discussion section that the NV center should be negatively charged. 

In the NV defect, three carbon atoms have sp 3 dangling bonds near the vacancy and 
three back bonds each to the lattice, while the nitrogen atom has also three back bonds 
and one dangling bond pointing to the vacant site. Since nitrogen has five valence electrons 
the negatively charged NV defect has altogether six electrons around the vacant site. The 
structure of the NV defect, including the definition of the symmetry <111> axis, is depicted 
in Fig. [H 
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FIG. 1: (Color online) The NV center viewed in perspective (left) and along the <111> direction 
(right); this direction is the C3 symmetry axis of the C3„ symmetry group of the defect. The vacant 
site is indicated by a small pink circle and the neighboring C and N atoms by grey and cyan balls, 
respectively. The first neighbor atoms around the vacant site have sp 3 dangling bonds pointing 
toward the vacant site, which are labeled (Xj {i = 1,4), as in the scheme used in the group theory 
analysis. 

The group theory analysis of the six electron model has been worked out previously for 
this defect.— We summarize the results using our notation and conventions: Since it is 
known that the carbon and nitrogen atoms relax outward from the vacancy,— 1 ^ we assume 
that the overlap between the dangling bonds (o"i_4) is negligible, that is, OiOj = 5ij. o~i_4 
are transformed under the operation of C% v point group forming the following orthonormal 
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FIG. 2: (Color online) The calculated spin-resolved single-electron levels with respect to the valence 
band maximum (VBM) in the ground state of the NV defect. Valence and conduction bands of 
the host crystal are shown as blue and pink shaded regions, respectively. The levels are labeled as 
in Eq. ((TJ) and their occupation is given for a negatively charged defect (a total of 6 electrons). 
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where < a < lis a parameter that determines the extent to which the nitrogen dangling 
bond is mixed in the 0i and <p2 defect states. There are two fully symmetric one-electron 
states (ai) and one doubly degenerate e state, with a total occupation of 6 electrons. We 
note here that the dangling bond of nitrogen is not mixed in the e state but only in the a\ 
states. It was found by Goss et a/.— using ab initio molecular cluster calculations that the 
two a\ states are lower in energy than the e state. As a consequence, four electrons occupy 
the di states and two electrons remain for the e state. Our calculated one-electron levels 
obtained by ab initio supercell calculations are shown in Fig. [2j As can be seen from this 
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analysis, the natural choice is to put the two remaining electrons in the e level forming an 
S=l state (by analogy to Hund's rule for the p-orbitals of the isolated group IV elements in 
the Periodic Table). 

In the C 3v point group the total wavefunction has 3 A 2 symmetry with S=l. In our 
special case we choose Ms = 1, so both electrons are spin-up electrons in the e level. 
The Csv symmetry can also be maintained by other occupations of the states. Putting 
two electrons into four possible quantum states of the degenerate e level, we end up with 
six possible multiplets (including the degeneracy): 3 A 2 , 1 A 1 , and 1 E. By taking the 3 A 2 
to be the ground state of the defect (as experiments indicate) there is no allowed optical 
transition to first order, since the spin state cannot be changed in a PL process. The 0i 
level is relatively deep in the valence band, so to a good approximation we can assume that 
it does not contribute to the excitation process. However, the <p 2 level in the gap is not very 
far from the e level. If one electron is excited from 2 into the e level (03 or 4 ) then either 
a 3 E or a l E multiplet is obtained. If both electrons are excited from <p 2 to 03 and 04 then 
a fully symmetric l A\ state is obtained. The only allowed transition is 3 A 2 ^ 3 E to first 
order. The electronic configurations of these states are explained in Table [H Most of these 



states were discussed in Ref 
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The Ms = ±1 triplet states can be described by a single Slater-determinant. However, 
the singlet states (except for the last 1 A 1 state) can be described by a linear combination 
of two Slater-determinants. The two singlet single-Slater-determinant states of the a^(2)e 2 
configuration are: 

|<T2<f2<7 3 <f 4 ) = ^[^(0, y) + 3 A 2 {0)) (2a) 

W 2 a 2 a 3 a 3 ) = ^E(0,x) + 1 ^(0)] (2b) 

where 1 E(0, A) is the multideterminant wavefunction of the singlet E state of Ms=0 at 
A-row (A = x or y) in Table [B The singlet single-Slater-determinant state of the a\(2)e 3 
configuration is: 

|<fa*3*4<f4) = ^[ 3 £(0, x) + 1 E{0, x)} (3) 

In summary, the following many-body states must be considered: 3 A 2 , 1 y4 1 , 1 E, 3 E, 1 E, 
and 1 v4i. The two triplet states are orthogonal to each other while the two l A\ and the two 
X E states theoretically can be mixed with each other. This will be discussed below. Our 
computational method described in the next section cannot take the spin-orbit and spin-spin 



TABLE I: The electronic configurations and the possible total wave functions with symmetry. 
For simplicity we abbreviate 01 — ► 2, etc. in the last column. Overbar in a wavefunction means 
spin-down electrons, the rest are spin-up electrons. We assume that ai(l) : o\ is fully occupied, so 
we do not show that part of the wave function here. In the second and third columns we give the 
symmetry of the total wavefunction (r) and its spin projection (Ms), respectively. In the case of 
doubly degenerate representations (E states) we designate which transforms as x or y in the last 
column. 
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interaction into account. From the energetic point of view those effects are marginal (within 
few meV) but they could have important consequences on the possible optical transitions 
and the spin state of the NV center.— 



7 



III. COMPUTATIONAL METHOD 



We use density functional theory with the local spin density approximation (DFT-LSDA) 
of Ceperley- Alder— as parameterized by Perdew and Zunger.— We employed three different 
codes and somewhat different methodologies to carry out the calculations. The geometry 
of the defect was optimized with the VASP codecs and the SIESTA code.^ The lat- 
ter utilizes numerical atomic orbitals with Troullier-Martins pseudopotentials. 41 We applied 
the high level double-^ plus polarization functions for both carbon and nitrogen atoms. In 
the SIESTA calculations no symmetry restriction was applied. The linear combination of 
atomic orbitals (LCAO) analysis of the defect states is straightforward in this methodology 
through the wavefunction coefficients obtained directly from the SIESTA calculations. In 
the VASP calculations we use a plane wave basis set with cutoff 420 eV(~30 Ry) which is 
adequate for well converged calculations using projected augmented wave (PAW) pseudopo- 
tentials for the C and N atoms.— ^ In the VASP calculations we applied the C^ v symmetry, 
and the energy of the ground state as well as that of the excited states are calculated by 
setting the appropriate occupation of the defect states in the gap. In the geometry optimiza- 
tion calculations, all the atoms were allowed to relax until the magnitude of the calculated 
forces was smaller than 0.01 eV/A. The hyperfine tensor of the NV center was calculated 
by the CPPAW code.— In the CPPAW calculations we used a 30 Ry cutoff for the plane 
wave basis with PAW projectors which is virtually equivalent with the methodology used in 
the VASP calculations. 

Convergence of calculated defect properties with supercell size is an important consider- 
ation. For this reason, we have chosen to model the NV center by using a large 512-atom 
simple cubic supercell. The lattice constant of the supercell (wl4.2 A) is four times larger 
than the lattice constant of the conventional cubic cell of diamond, ao=3.54 A. We used the 
T-point sampling in the Brillouin-zone which corresponds to sampling finer than a 6 x 6 x 6 
grid of the primitive lattice; this provides a well converged charge density. It is also ad- 
vantageous to restrict the calculations to the T-point in order to keep the degeneracy of 
the e defect states which may split in a general k-point sampling of the Brillouin-zone. We 
checked that the geometry was practically identical (to within 0.01 A) going from a 216- 
atom fee supercell to the 512-atom simple cubic supercell. We will show that the calculated 
spin density, for instance, decays at much shorter distance than the lattice constant of the 
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supercell. Thus, the 512-atom supercell is adequate to represent the isolated NV defect in 
a realistic manner. 

We calculated the hyperfme tensor of the defect with the optimized geometry obtained 
by the VASP code. In the calculation of the hyperfme tensor the relativistic effects are 
taken into account.— The hyperfme tensor of nucleus I consists of the Fermi-contact term 
(first parenthesis in the following equation) and the dipole-dipole term (second parenthesis): 

4> = ±f d»r n,(r hj7eR *[(f f W ) + - |)] (4) 

where n s (r) is the spin density of the spin state S, 7/ is the nuclear Bohr-magneton of 
nucleus / and 7 e is the electron Bohr-magneton. The Fermi-contact term is proportional to 
the spin density localized at the place of the nucleus which is dominant compared to the 
dipole-dipole term. The ratio of the Fermi-contact and dipole-dipole terms characterizes the 
shape of the spin density. 

The contribution of s-like wave functions to the charge density has a large effect on 
the Fermi-contact term but negligible effect on the dipole-dipole term, since the s-like wave 
function has a maximum at the positions of the nuclei and it is an even function. In contrast 
to this, the contribution of p-like wave functions to the charge density has a negligible effect 
on the Fermi-contact term but a large effect on the dipole-dipole term, since the p-like 
wave function has a node at the place of nuclei and it is an odd function. Typically, the 
contribution of the dipole-dipole term is significant for the spin density built from well- 
localized dangling bonds, that is, the sp 3 hybrid orbitals (see Table HB . We note that the 
pseudopotential methodology produces artificially smooth wave functions close to the nuclei, 
therefore only the all-electron PAW methodology can provide reliable hyperfme tensors. In 
the PAW methodology the calculation of the hyperfme tensor is somewhat more subtle than 
Eq. PQ) shows and the dipole-dipole term is not fully calculated (see the appropriate note in 
Ref J45l) which causes about 0.3 MHz inaccuracy in the calculated dipole-dipole term in our 
case. This error could be important to take into account for hyperfme tensors with small 
matrix elements (<3 MHz). The total spin density, n s (r), can be defined as 

n s (r) = n up (r) - n down (r) 

where n up (r) and n down (r) are the spin densities built from spin-up and spin-down electrons, 
respectively. Taking the M s =l state of the 3 A 2 state, as shown in Fig. (2J we expect that 



n s (r) will be positive. We will show that this is not true for the entire space around the 
defect. 

IV. ATOMIC AND ELECTRONIC STRUCTURE OF THE NV DEFECT 
A. Geometry and electronic levels 

We begin with a discussion of the geometry optimization of the negatively charged NV 
center obtained from the SIESTA calculations using spin-polarization and no symmetry 
restrictions. The defect automatically finds t the S=l state and maintains the symmetry 
of the original, unrelaxed structure. The calculated one-electron defect levels are shown in 
Fig. [21 The first neighbor C and N atoms clearly relaxed outward from the vacancy. The 
calculated distances from the vacant site are 1.63 and 1.69 A for the C atoms and the N 
atom, respectively, that is, the N atom relaxes more than the C atoms. We note that the 
C-vacancy distances are the same within 0.0002 A without any symmetry constraints after 
geometry optimization. The N-C bond lengths are 1.46 A, while the bond lengths of C 
radicals are 1.50 A, which is not far from 1.44 and 1.45 A, respectively, obtained in an LDA 
molecular cluster calculation.— The localized basis sets can have problems in the description 
of vacancy-like defects. Our VASP calculation, which employs a plane wave basis set, shows 
basically the same geometry after optimization: the calculated distances from the vacant 
site are 1.62 and 1.68 A for the C atoms and the N atom, respectively. Thus, we conclude 
that the double-^ plus polarization basis provides results very close to those of the converged 
plane wave basis set. 

In an earlier work by by Luszczek et al.,— a plane wave basis set with pseudopotentials 
was employed in a 64-atom cubic supercell using 2x2x2 Monkhorst-Pack Brillouin Zone 
sampling 46 to investigate the NV defect in diamond. In that work, only the nearest neighbor 
atoms to the vacant site were allowed to relax without symmetry restrictions and a geometry 
close to C% v symmetry was obtined; the largest deviation in the C-vacancy distances was 
about 0.001 A. The calculated distances from the vacancy were 1.67 and 1.66 A for the C 
atoms and the N atom, respectively, which shows the opposite trend from what we find both 
in the SIESTA and in the VASP calculations in the larger unit cell. Most probably the 
difference is due to the insufficient relaxation condition restricted only to the first neighbor 
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atoms around the vacant site. 

We plot the wave functions of the defect states obtained by the LSDA calculations in 
Fig. [3J The group theory analysis based on the single-particle picture describes very well 

<m> ~r (in) 
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FIG. 3: (Color online) Isosurfaces of the calculated wavefunctions of the ai(l), ai(2), e x , e y defect 
states, shown in side (left) and top (right) views relative to the <111> axis. Blue (red) isosurfaces 
correspond to negative (positive) values of the wavefunction. The small pink circle represents the 
position of the vacant site, while the grey and cyan balls show the C and N atoms, respectively. 
We show the atoms up to the second neighbor from the vacant site. 



the defect states. Naturally, the defect states are not strictly localized on the first neighbor 
atoms of the vacancy but the largest portion of the wave functions can be indeed found 
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there. The SIESTA calculation gives a ~ 0.7 for the parameter that appears in Eq. ([I]). 
This means that the N orbital is mostly localized on the ai(2) : <p 2 defect level and has no 
amplitude on the e levels. Therefore, the nitrogen atom is only very weakly spin-polarized in 
the 3 A 2 state, while it is strongly spin-polarized in the 3 E state (when one electron is excited 
from the ai(2) to level to the e level). This is clearly shown on Fig. HI It is apparent, from 
this figure, that the N atom is only weakly polarized (small negative spin density) in the 
3 A 2 state while it is strongly polarized in the 3 E state comparable to the C ligands (large 
positive spin density). The spin density is always highly localized on the three C atoms 
around the vacant site (orange lobes in the figure). 




FIG. 4: (Color online) Calculated spin density isosurfaces in the Ms=l state for the 3 A 2 (left) and 
the 3 E state (right). The vacant site is depicted by a small pink sphere at the center of each plot. 

We also checked the situation when we optimized the geometry with the condition S=0. 
We already showed that the \<j 2 <j 2 <J3<J3) state is not an eigenstate with symmetry. In 
addition, this state is a Jahn- Teller unstable system. Indeed, the defect reconstructs to C\h 
symmetry to remove the degenerate e level. However, this configuration is about 0.3 eV 
higher in energy than the 3 A 2 state with C% v symmetry. 



12 



B. Energy sequence of multiplets 



From the structural analysis, we conclude that the dangling bonds around the vacancy 
do not form long bonds which could be the driving force of the reconstruction. Instead, 
the atoms relax outward from the vacancy and retain the strongly localized dangling bonds 
pointing to the vacant site which maintains the C3,, symmetry. Since the degenerate e defect 
level is only partially occupied, this is a typical situation where configurational interaction 
plays a crucial role. As shown above, most of the singlet eigenstates can be described only 
by multi-determinantal wavefunctions in C^ v symmetry. The optical transition takes place 
between the triplet states. We already showed the results on the 3 A 2 state. In the VASP 
calculation it is possible to set the occupation of one-electron states. The 3 E state can be 
achieved by setting zero occupation for the spin-down 02 level and full occupation of 03 
spin-up and spin-down levels. The energy of the 3 E state can be calculated in the fixed 
geometry of the 3 A 2 state which yields the vertical ionization energy. Upon the excitation 
of the electron the nuclei can relax to find the minimum energy in the new configuration 
space. This relaxation can take place with the help of phonons around the defect. The ZPL 
transition corresponds to that energy where phonons do not participate, between the energy 
minima of the two configurations as shown in Fig. [5l 




FIG. 5: Energetics of photoluminescence absorption: VI is the vertical ionization energy, ZPL is 
the zero-phonon line transition, E is the total energy and Q the configuration coordinate. 

The calculated vertical ionization energy is 1.91 eV within LSDA. We found that the NV 
defect significantly relaxes due to this internal ionization. The C- vacancy distance is 1.67 A 
while the N-vacancy distance is 1.61 A in the 3 E state. This shows the opposite trend than 
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what was found in the 3 ^2 state. This may be understood as follows: the N atom is strongly 
spin-polarized in the 3 E state compared to 3 A2 state, while the C ligands will be somewhat 
less spin-polarized which induces different charge transfer between the atoms in the 3 E state, 
and leads to a different geometry. The calculated relaxation energy (the Fracnck- Condon 
shift) is 0.2 eV. From this, we find a ZPL energy of 1.71 eV which can be tentatively 
compared to the experimental value of 1.945 eV.— We note that a similar value (1.77 eV) 
was found by the LSDA molecular cluster calculation. 34 The LSDA excitation energy and 
the experimental transition energy are remarkably close to each other. This shows that the 
self-interaction error of LSDA for these defect levels does not differ too much, which is not 
unexpected since both of the defect states are basically valence band derived (from sp 3 -like 
hybrid orbitals). Nevertheless, the calculated Franck-Condon shift, which is the relaxation 
energy defined as the energy difference between the vertical ionization energy and the ZPL 
energy, should be even more accurate than the calculated internal ionization energy. Indeed, 
the PL spectrum shows a broad phonon spectrum even at low temperature, and the intensity 
of the ZPL line is relatively small compared to the phonon side bands which indicates a large 
Franck-Condon shift. 

In addition to the triplet states it is worthwhile to calculate the energies of the singlet 
states because they play a significant role in the emission process but these states have not 
yet been measured directly in experiments (see Refj23l and references therein). LSDA is not 
a suitable methodology to calculate these energies accurately. Beside the self-interaction 
error (which is relatively small for these defect levels as discussed above), LSDA gives the 
charge density of the interacting electrons, which is expressed in terms of the non-interacting 
Kohn-Sham particles. In the Hartree-Fock language the total wave function, which is not 
directly used in DFT-LSDA, is expressed by a single Slater-determinant. Thus, the usual 
DFT-LSDA calculation cannot represent most of the singlet eigenstates of the NV center. 
This can be corrected by using a perturbation theory within DFT-LSDA as explained by 
Lannoo et al.. Ar Instead, we adopt here a more approximate method 4 ^ which was already 



applied for some states of this defect in RefyJ. Our goal is to give approximate energies of 



the singlet eigenstates in order to predict their sequence, and we do not attempt to provide 
energies directly comparable to the experiments. 

One of the 1 A 1 states can be described by a single Slater-determinant (when the ai(2) level 
is totally empty and the e level is fully occupied). This state can be calculated with geometry 
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optimization by LSDA as explained for the 3 E state. We therefore concentrate on the 
remaining 1 Ai and two 1 E states. While the \a 2 a 2 a 3 a 3 ) , l^^&s&i) , and \cr2cr3cf3cf4) single 
Slater-determinants are not eigenstates, they can be expressed as the linear combination 
of different eigenstates shown in Eqs. (|2ap . (pbj) and ([3]). von Barth^ has shown that the 
LSDA total energy {£ ) of the mixed state can be expressed as the appropriate sum of the 
energy of the eigenstates, which yields the following equations: 

£[\a 2 a 2 a 3 a 4 )] = \(S?E\ + £[ 3 A 2 ]) (5a) 

e[\W2<rz<f z )\ = \{BM + £[ l A x ]) (5b) 

£[\a 2 a 3 a A d 4 )\ = ^(£[ 3 E] + £[ l E]) (5c) 

In Eqs. ( )5a|l and ([3c]) we assume that the energy of the triplet states with Mg = and Mg = 1 
is the same. This is a very good approximation since, for instance, the experimentally 
measured splitting is about 2.88 GHz (few /xeV) for the ground state due to spin-spin 
interaction,— while the spin-orbit splitting for the 3 E state is expected to be within few 
meV, which is far beyond the accuracy of LSDA calculations. The energy of the mixed states 
on the left hand side of Eqs. (!5al) -(l5cl) can be calculated directly by LSDA. Since £[ 3 ^b] is 
known, ^f 1 ^] can be determined from Eq. fl5a|) . Similarly, ^f 1 -^] can be determined from 
Eq. (I5cj) . By combining Eqs. (I5a|) and (15bl) we arrive at: 

E[ l A x ] = 2(£[\a 2 a 2 a 3 a 3 )} - £[\a 2 a 2 a 3 a 4 )}) + £[ 3 A 2 }) (6) 

Thus, the l A x state can be also determined. 

von Bartb— applied this approach successfully to calculate the energy of atoms in different 
states. Formally, this method can also be applied to the NV defect in diamond, but attention 
must be paid to relaxation effects. Generally, if the electron state changes then it may 
imply also relaxation of the ionic positions, as was the case for the 3 E state discussed 
earlier. Relaxation effects cannot be taken into account with this methodology, since only 
the energy of the mixed states can be calculated directly by LSDA, and the relaxation 
of the mixed state is meaningless. In other words, the geometry must be fixed in these 
calculations. The occupation of the e state varies in the case of the l A x and 1 E states of the 
a\(2)e 2 configuration. We assume that the geometry would involve negligible change from 
the geometry of the 3 A 2 state which belongs also to the af(2)e 2 configuration. Therefore, 
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we fix the geometry obtained in the 3 A 2 state in the process of calculating the x Ai and 1 E 
states of the a\(2)e 2 configuration. Using the same argument we fix the geometry obtained 
for the 3 E state in the calculation of the l E state of the a\(2)e 3 configuration. With these, 
we find the following energy sequence of the multiplets: 

si \ i \ a -I ~0.0 eV /-r1 i n ~0.9 eV ,-.ri m ~0.8 eV nrHi-n ~0.5 eV ,-.1-1 ,-n asl.3 eV /->r1 i n 

£[ A 2 ] > £[ Ax] > £[ E] > £[ E] > £[ E] > £[ A^, 

that is, the deeper 1 Ai state is close in energy to the 3 A 2 state, and the deeper 1 E state is 
below 3 E state. The energy differences between the triplet states were already discussed. 

We focus next on the singlet states. Within our approximate methodology the 3 A 2 
and l A\ states are almost degenerate. Close inspection of Eq. ([6]) reveals that the energy 
sequence of the 3 A 2 and 1 Ai states depends on the energy difference of two singlet states: 
{£\\cr 2 a 2 a^(7z)] — £[\ct 2 <j 2 <j 3 (J4)]). We obtain almost zero for this energy difference. However, 
the LSDA self-interaction error may be larger for the \cr 2 <f 2 a3(f3) state (where the exchange 
energy of 0-3 appears) than for the \<j 2 <j 2 <j state (where the electrons occupy spatially 
orthogonal orbitals). This may raise the energy of the ^P^i] state. Nevertheless, this can 
be partially compensated by the relaxation effect of the 1 A 1 state which we are neglecting 
by necessity as explained before. An additional issue is the possible mixing with the higher 
1 Ai state. We argue that the two X A\ states are not likely to mix because they are very far 
from each other in energy. 

We conclude from the above analysis that the 1 A 1 state is indeed close in energy to the 
3 A 2 state. This may imply a very complicated fine structure of the states. Taking the spin- 
spin interaction into account the 3 A 2 state splits to A^Mg = 0) and E(M$ = ±1) states, 
while 1 Ai becomes Ax (Ms = 0) (the spin state is not a good quantum number anymore, 
just its projection). If the two Ai (Ms = 0) states are close in energy then they may mix 
with each other. It is also important to notice that the energy of the deeper l E state falls 
between the energies of the triplet states. The relaxation effect may lower the energy of this 
state. In addition the two 1 E states are not far from each other in energy, so the off-diagonal 
elements in the hamiltonian may not be neglected. The mixing of the two X E states would 
further lower the energy of the deeper l E state and would raise the energy of the higher 1 E 
state. The final conclusion is that there are two singlet states between the triplet states, and 
the 1 A\ state is much closer in energy to the ground state than the X E state. The energies 
of the other two singlet states are certainly above that of the 3 E state. 
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The deeper singlet l A\ and 1 E states can play an important role in the emission process 
of the NV center. Experiments indicate that there should be a possibly long-living singlet 



state between the triplet states (see Ref ]23l and references therein). Usually, the singlet 1 Ai 
is considered in this process. However, our calculations indicate that there are two singlet 
states between the triplet states. Goss et a/. 34 reported the sequence of 3 A 2 , 1 E, 1 Ai, and 3 E 
states which is surprising in light of the previous discussion, as these states were obtained by 
LSDA molecular cluster calculations. Manson et alM recently showed that if the 1 E state 
is above the 1 A\ state, the known properties of the emission can be consistently explained 
as with the original singlet 1 A\ model, with the only difference being that the l E and l A\ 
states both contribute to the spin polarization process during the optical cycling and that the 
effect will be more efficient. Having the l Ai state higher in energy than the l E state would 
result in no change in spin orientation during optical cycling which is in contradiction with 
experiment.— From this point of view our results obtained from approximate calculations 
are consistent with the PL experiments. 

V. HYPERFINE CONSTANTS 

As mentioned in the Introduction, the NV center in diamond is a promising candidate 
to realize qubit solid state devices operating at room temperature (RefJl4j and references 
therein). The qubit is the non-zero (S=l) electron spin ground state which can interact 
with the neighbor 13 C isotopes possessing 1=1/2 nuclear spin via hyperfine interaction. The 
natural abundance of the 13 C isotope is about 1.1%, so we can assume the same abundance 
in the diamond lattice. In conventional EPR measurements, the EPR absorption signal 
is detected on the ensemble of the defects in the diamond sample. The sample should be 
thick enough for absorption measurements and the concentration of the defects should be 
sufficiently high. Finally, a large number of defects is measured at the same time by EPR, 
so statistics can be applied to analyze the data. If the spin density is strongly localized on 
three symmetrically equivalent C ligands (see Fig. H]) then the probability of finding one 13 C 
atom among them is given by the binomial distribution and is about 3.2%. Finding two 
or three of them has negligible probability. Due to the 1=1/2 nuclear spin, the hyperfine 
interaction splits to two lines with I z = 1/2 and I z = —1/2 and therefore the intensity ratio 
between the main hyperfine line (involving no hyperfine interaction with 13 C atoms) and 
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the hyperfine line associated with the C ligands will be roughly 1.5%. This makes the EPR 
measurement on 13 C hyperfine interaction a challenging task, since the signal to noise ratio 
should be very good and the intensity of the EPR signal should be strong enough and stable 
to identify the satellite hyperfine lines due to 13 C isotopes. 

The 13 C hyperfine interaction has been detected in EPR from two sets of three symmetri- 
cally equivalent C atoms by Loubser and van Wyk.— ^ The larger hyperfine constants were 
associated with the C-ligands of the NV center. In addition, the hyperfine interaction of 
the M N isotope was found in NV center by EPR and ENDOR measurements.— 1 ^ To our 
knowledge, hyperfine interaction with other 13 C isotopes has not been measured directly by 
EPR. We can estimate the localization of the charge density on the C and N atoms that are 
immediate neighbors of the vacant site, from the linear combination of the atomic orbitals 
as they appear in the wavef unctions. He et a/.— estimated that 72% and 0.2% of the charge 
density is localized on the three C ligands and the N atom, which leaves approximately 28% 
of the charge density to be spread in the lattice. Wrachtrup et a/. 3 speculated that the spin 
density decays exponentially as a function of the distance from the vacant site. Since the 
hyperfine constants are roughly proportional to the spin density, as Eq. ((H])) indicates, these 
authors proposed that 9 or more carbon nuclei should have a hyperfine value of 70 MHz 
in the second neighborhood while the more distant carbon atoms should have a hyperfine 
constant smaller than 10 MHz. This proposal is not entirely consistent with the known 
EPR data, since isotropic hyperfine splitting of 5.4 G(~15 MHz) was measured from three 
13 C isotopes 4 - whereas a value of 70 MHz hyperfine splitting should be measurable by EPR 
because it would not be obscured by the main EPR line. 

In their theoretical treatment of the NV center, Luszczek et alr^- claimed that they can 
support the proposal of Wrachtrup et al— based on their ab initio results. They optimized 
the geometry without symmetry constraints only for the first neighbor atoms of the vacancy 
in a 64-atom supercell. The adequacy of this restriction was already discussed above; with 
this restriction, while the C^ v symmetry is almost retained, the calculated Fermi-contact 
hyperfine interactions for the three C ligands deviate from each other more than 10% (see 



Table 2 in Ref. 
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3) . This suggests that the spin density was not adequately converged in that 
calculation since the small deviation in the geometry from the symmetry could not imply 
such a large discrepancy in the calculated hyperfine field. These authors also calculated the 
Fermi-contact hyperfine interaction for the C atoms situated about 2.5 A away from the 
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vacant site. The reported numbers were about an order of magnitude smaller than for the C 
ligands, which led to their conclusion that the spin density and the corresponding hyperfme 
constants should decay fast for other C-atoms farther from the vacant site. 

In addition to the problem of the inconsistent values of the hyperfme interaction for the 
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must be mentioned: 



three C-ligands, several other issues related to the results of Ref . 

(i) These authors reported actually the hyperfme field and not the hyperfme constant. How- 
ever, the conversion from the hyperfme field, which is the magnetization density on the atom 
and it is a number directly obtained from the computation, to the hyperfme constant is not 
unique; therefore, it is very difficult to compare the calculated values to the experimental 
data. 

(ii) Only the Fermi-contact term was calculated while the dipole-dipole term can be also sig- 
nificant; this is known to be the case for the C ligands from experimental measurements.—^ 

(iii) The hyperfme interaction with distant C atoms could be very important for qubit appli- 
cations based on this defect, so the hyperfme tensor must be calculated at larger distances 
from the vacancy. 

Based on these arguments, we believe that the nature of the spin density and the corre- 
sponding hyperfme interaction with the 13 C isotopes has not yet been explored in detail 
despite of its high importance. 

Toward establishing the nature of this interaction, we show first the calculated spin 
density in our 512-atom supercell in Fig. [6j As expected, the spin density is highly localized 
around the three C-ligands nearest to the vacant site (orange lobes in the figure), within a 
radius of 1 x ao from this site. The spin density practically vanishes at distances > 2 x ao 
from the vacant site. Below the three C atoms (see Fig. [6]) there is the N atom with a small 
negative spin density. Some C atoms farther from the vacant site also have negative spin 
density. The spin density extends mostly on a plane perpendicular to the (111) direction and 
no measurable spin density can be found below the N atom. On the N atom, the spin density 
is negative. It was shown earlier that the spin density comes mostly from the spin-polarized 
e level localized on the C dangling bonds, and due to symmetry reasons, orbitals related to 
the N do not appear in the e level. The tiny negative charge density on the N atom can be 
explained by the polarization of its core states: since the nuclear Bohr-magneton of the N 
atom is positive, the Fermi-contact term will be negative (see Table HI) . Overall, the spin 
density is spread on a plane perpendicular to the (111) direction. There are some C atoms 
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FIG. 6: (Color online) (a) The 512-atom cubic supercell with all C atoms shown (atoms within a 
64-atom supercell are shown in red). The suprercell size in A and in units of the conventional cubic 
cell lattice constant, do, is indicated on the side of the cube, (b) Perspective view of the calculated 
spin density isosurfaces in the Ms=l state. Only atoms up to the second nearest neighbor of the 
vacant site are shown, (c) View along the (111) direction indicating the C^ v symmetry of the 
spin density, which is given in colored contours. The black lines denote the size of the 512-atom 
supercell, while the red dotted ones show the boundaries of the 64-atom supercell. 

which have significant negative spin density which results in negative Fermi-contact term. 
No measurable hyperfine interaction (spin density) can be found for the C atoms below 
the N atom. Loubser and van Wyk^ speculated that the 15 MHz 13 C isotropic hyperfine 
splitting comes from the three C atoms bonded to the N atom. Our calculation negates this 
possibility. 

Table HI1 shows that there is a local maximum of the spin density at R vac ~3.9 and 5.0 A, 
where R vac is the distance from the vacant site. At R vac >6.3 A the calculated hyperfine 
constants are below 1 MHz which means that the spin density vanishes at R vac ~ 2xa . 
Apparently, the 64-atom supercell is too small to capture these properties due to the artificial 
overlap of the spin density caused by the periodic boundary conditions. 

We also find that the charge density does not decay monotonically like an exponential 
function, but behaves more like a wavepacket, that is, as f(x) = sm(x)/x (see Fig. EI) . 

It is important to compare the calculated hyperfine values with the known experimental 
data in order to estimate the accuracy of our calculations. For that purpose, the hyperfine 
constant of the N atom is not the best choice because its value is very small and it is 
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TABLE II: The calculated principal values of the hyperfine tensor (columns 3 to 5) compared 
to the known experimental data (columns 6 to 8) in MHz. The average of the three principal 
values yields the Fermi-contact term. The difference between the principal values and the Fermi- 
contact term gives the dipole-dipole term. Only atoms with a signal larger than 2 MHz are shown. 
The symmetrically equivalent number of C atoms is shown in the first column and their distance 
in A from the vacant site in the second column. The experimental data on 14 N is taken from 
Refs. 
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49l . EPR studies can directly 



J49j. Experimental data on 13 C atoms were taken from Ref ■ 
measure only the absolute value of the hyperfine constants, which is indicated by adding a ± sign to 
experimental values. The calculated hyperfine tensors can be used for comparison with spin-echo 
measurements (see text). 



Atom 


Rvac 


A n 


A 22 


A33 


A ex V 


A ex 'P 
A 22 


A exp 
A 33 


14 N 


1.68 


-1.7 


-1.7 


-1.7 


±2.1 


±2.1 


±2.3 


3 C 


1.61 


109.5 


110.2 


185.4 


±123 


±123 


±205 


6 C 


2.47 


-4.8 


-3.7 


-1.5 








3 C 


2.49 


-7.4 


-7.3 


-5.8 








6 C 


2.90 


2.8 


3.3 


4.6 








3 C 


2.92 


1.4 


2.4 


2.9 








3 C 


2.93 


3.4 


4.7 


4.9 








6 C 


3.85 


13.5 


14.2 


19.4 








3 C 


3.86 


12.8 


12.8 


18.0 


±15.0 


±15.0 


±15.0 


6 C 


4.99 


2.6 


2.7 


3.8 








3 C 


5.00 


1.5 


1.5 


2.2 









caused by only indirect spin-polarization effect. As mentioned above, there is an inherent 
inaccuracy in the calculated dipole-dipole term of about 0.3 MHz, therefore it is reasonable 
to consider only values of the hyperfine constants that are significantly larger than this 
limit. Accordingly, we restrict the comparison to values that are higher than 2 MHz. By 
comparing the hyperfine constants of the C-ligands we estimate the inaccuracy for both the 
Fermi-contact term and the dipole-dipole term to be about 10%. This is usual for LSDA 
calculations.— JJL22i2i The calculated ratio between the Fermi-contact term and the dipole- 
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FIG. 7: Variation of the principal values of the hyperfine tensor An and A33 as a function of the 
distance from the vacant site (R vac ) for the set of 3 symmetrically equivalent C atoms (3C). In the 
inset we show changes farther from the vacant site on a finer scale. The variation for the set of 6 
symmetrically equivalent C atoms follows closely the one shown for 3C. 

dipole term agrees remarkably well with experiment for the C-ligands. This shows that 
the shape of the spin density is very well described by LSDA. This ratio indicates that the 
p-functions dominate by about 90% in the dangling bonds, so those are more p-like orbitals 
than sp 3 hybrids. The reason is most likely the outward relaxation of the C-atoms from 
the vacant site. The plotted wavefunctions indeed show a very strong p-contribution of the 
e levels in Fig. [3] while the spherical s-contribution is very small. The shape of the wave 
functions can be directly compared to the proposed wave functions in Eq. ([I]) from the group 
theory analysis of the defect diagram. 

An additional 15 MHz 13 C hyperfine splitting was measured by Loubser et alA- using EPR 
but the measured spectrum of the NV center was not shown in detail, therefore we cannot 
comment on the accuracy of this measurement. Nevertheless, that work stated that the 
hyperfine splitting is isotropic, that is, the dipole-dipole term is negligible, and the relative 
intensity of the hyperfine satellite line and the main EPR line indicates the involvement of 
three symmetrically equivalent C atoms.— The likely candidate for this signal is found at 
the third neighbors of the vacant site at R vac =3.86 A. The spin density of these atoms is 
shown as yellow lobes above the small green lobes at the highest position on the side view 
in Fig. [HI The calculated anisotropy of this hyperfine interaction is about 3-4 MHz. If the 
uncertainty in the measurement is in this range due to line broadening, then the signal could 
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be detected as nearly isotropic. However, our calculations reveal that six additional C atoms 
have similar hyperfine splitting at the third neighbor distance of R vac =3.85 A, corresponding 
to the six yellow lobes laterally spread farthest from the vacant site, which is most obvious 
from the view along the (111) direction in Fig. El This means that the hyperfine splitting 
due to these six C atoms would have to be detected simultaneously with the other three 
ones. Since the difference in the hyperfine splitting of these two sets of atoms is small, the 
six-atom set could obscure the signal of the three-atom set showing an effective relative 
intensity associated with six symmetrically equivalent C atoms. We suggest that this part 
of the spectrum should be re-investigated in detail experimentally. The hyperfine splitting 
of 7-8 MHz from the 13 C atom may be also detectable by EPR while the other signals may 
be too small and hence obscured by the main hyperfine lines. 

Recently, individual NV centers have been detected by spin-echo measurements.— 1 ^ In 
particular, detailed results for six NV centers in diamond have been reported.— The spin-echo 
measurements have detected the coherent state of the electron spin coupled with a proximal 
13 C nucleus. The coupling is due to hyperfine interaction between the electron spin and 
the nuclear spin of 13 C isotopes. The resulting spin-echo signals show a rapidly oscillating 
function enveloped by a more slowly oscillating function.— These authors proposed a theory 
to explain this signal, and they concluded that the fast modulation frequency is due to the 
effective magnetization density of the electron spin felt by the 13 C nucleus, which is the same 
as the hyperfine interaction. The modulation frequency can be well approximated as the 
norm of the hyperfine tensor projected to the symmetry axis, which leads to the following 
expression within our formulation of the problem: 



= V([ui]*4f) 2 ( 7 ) 

where [111] is the appropriate projection vector and A\- is the hyperfine tensor of nucleus 
X. Since we calculate the full tensor (see Eq. PJ), the calculation of the modulation 
frequency is straightforward. Modulation frequencies have been reported for four single NV 
defects (see Fig. 4B in RefJl3l). A single modulation frequency was measured for each NV 
center at ^2 MHz, ^4 MHz, ^9 MHz, and «14 MHz, respectively, so these particular NV 
centers had hyperfine interaction with one 13 C isotope in the lattice.— It would be useful 
to compare the hyperfine interaction measured by EPR and spin-echo techniques which can 
be an additional validation of the theory developed by Childress and co-workers.— As was 
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explained earlier, conventional EPR tools have limitations in detecting 13 C isotopes. A 13 C 
enriched sample would be useful to study this defect experimentally in more detail. 

The largest 13 C hyperfme splitting corresponds to a modulation frequency that is too 
large to be detected by spin-echo measurements. However, an isotropic 13 C hyperfme split- 
ting of 15 MHz has also been reported.— The isotropic signal means that the modulation 
frequency should be also about 15 MHz. This is very close to one of the measured modula- 
tion frequencies at ~14 MHz.— From the calculated hyperfme tensors in Table HT| this signal 
must originate from a C-atom which is a third neighbor of the vacant site at R vac =3.86 A. 
The calculated modulation frequency is ~16 MHz which is close to the measured one taking 
into account the limitations of our computational method. The 9 MHz spin-echo signal 
can originate only from the atoms at R vac =2.49 A, and contributions from other neighbors 
can be safely excluded. In this way, the origin of the signal can be identified. The 4 MHz 
spin-echo signal can originate either from atoms at R vac =2.90 A or at R vac =2.93 A taking 
into account computational limitations. In the first case, six symmetrically equivalent C 
atoms are involved, while in the second case a set of three symmetrically equivalent C atoms 
are involved. It is difficult, if not impossible, to identify the origin of the 1 MHz spin-echo 
signal, which is beyond the accuracy of our calculations. Nevertheless, the calculations in- 
dicate that this signal could arise from at least 12 C atoms. Most of them are far from the 
vacant site but some are closer than the atoms that give rise to the ~14 MHz signal, as is 
evident from Figure [3 

The number of symmetrically equivalent atoms is also specific to the individual hyperfme 
constant, and therefore to the modulation frequency. Because of the symmetry, sets of 
three or sets of six C atoms are equivalent with each other. The relative occurrence of the 
modulation frequencies measured by spin-echo experiments helps in identifying the equiv- 
alent atoms around the vacancy. Four samples represent a rather limited set of values for 
statistical analysis, so the relative occurrence of 13 C isotopes picked up by these measure- 
ments cannot be used for such analysis. A much larger number of NV samples is needed in 
the spin-echo measurements in order to use the relative occurrence of the 13 C isotope signals 
for the identification of individual atoms in the diamond lattice. 
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VI. SUMMARY AND CONCLUSIONS 



In this work we investigated the negatively charged nitrogen-vacancy center in diamond 
in detail by ab initio supercell calculations using density functional theory methods. We 
showed that the energy sequence of multiplet states is 3 A 2 , 1 A i , 1 E, 3 E, 1 E, and 1 A 1 . This 
means that the singlet 1 E state enhances the spin polarization process during the optical 
cycling of the defect. The center has non-zero spin ground state. The full hyperfine tensor for 
a large number of atoms around the defect was calculated for the first time. The calculated 
hyperfine constants of the C ligands agree well with the experimental values detected by 
electron paramagnetic resonance tools. However, there is a controversy about the number 
of symmetrically equivalent carbon atoms of the second highest hyperfine interactions when 
these are compared to experiment. We propose that that part of the electron paramagnetic 
resonance spectrum should be re-investigated in detail in order to clarify this issue. Our 
calculations reveal that the spin density of the ground state is spread in the lattice, mostly 
on a plane perpendicular to the (111) direction defined by the positions of the N atom 
and the vacant site, and that it does not decay monotonically from the vacant site. As a 
consequence, only a certain number of 13 C isotopes can interact with the electron spin, which 
can be used as qubits for quantum computing. Using the limited number of measurements 
that have been recently published for single nitrogen-vacancy centers detected by spin-echo 
measurements, we were able to identify some individual atoms around the defect. Our results 
contribute to the understanding of the spin-echo signals of the defect, which is a crucial step 
towards realization of the qubit concept in this system. Additional spin-echo measurements 
in NV samples will help identify other individual 13 C atoms around the defect. 
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